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ABSTRACT 

The efficiency and uniqueness of the diffusive shock acceleration is studied 
on the basis of the novel kinetic solutions. These solutions obtained earlier 
(paper I) selfconsistently describe a strong coupling of cosmic rays with the 
gas flow. They show that the dependence of the acceleration efficiency upon 
physical parameters is critical in nature. In this paper we investigate a steady 
acceleration in the parameter space formed by the injection rate is, the upper 
cut-off momentum p\ and the Mach number M while the flow compression R 
serves as an order parameter. We determine a manifold of all possible solutions 
in this parameter space. 

To elucidate the differences between the present kinetic results and the well 
known two-fluid predictions we particularly focus on the v — > 0, M — > oo limit 
where the two-fluid model suffers from especially serious closure problems and 
displays an 'unphysical' behavior. We show that in contrast to the two-fluid 
model three different solutions occurs also for arbitrarily large M provided 
that pi is sufficiently high. The three solutions appear together only if the 
injection rate v lies between two critical values, v\ < v < z/ 2 . For v < ui(M,pi) 
only the inefficient solution is possible. For v > u 2 (M,pi) only the efficient 
solution with a very high cosmic ray production rate occurs. On the basis of 
the obtained bifurcation surface R(v,pi) we consider the limit p\ — > oo, v — > 
which completely uncovers the long debating anomalies of the two-fluid model. 

The constructed steady state manifold that, at least partially is an attractor 
of a time dependent system, allows us to speculate on the nonstationary 
acceleration. 



Subject headings: acceleration of particles, cosmic rays, diffusion, hydrodynamics, 
shock waves, supernova remnants 
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1. Introduction 

The question, how efficient the diffusive shock acceleration may be, arose naturally 
when the first test particle calculations of this process became available (Krimsky 1977; 
Axford, Leer & Skadron 1977; Bell 1978; Blandford & Ostriker 1978; see |Drury 1983| ; 
Blandford fc Eichler 1987| and |Jones fc Ellison 1991| for a review). This is because the 
backreaction of accelerated particles (cosmic rays (CRs) in the astrophysical context) on the 
shock structure is very strong and leads normally to a significant increase of the compression 
ratio. Such an accelerating shock should therefore be thought of as a strongly nonlinear 
dynamical system with a pronounced selforganization. Neither the particle spectrum nor 
the hydrodynamic flow structure can be calculated independently. Furthermore, since the 
diffusion length of particles increases with momentum, particles with higher momenta 
sample longer parts of the shock transition. This makes any kind of moment description 
very difficult. However, the first essentially nonlinear calculations of this acceleration 
process were performed within the hydrodynamic approach. 



1.1. Success and limitations of the two-fluid model 

The above arguments suggest that the problem is kinetic in nature, which makes the 
usage of any fluid theory for describing the acceleration process questionable. At the same 
time quite a deep insight can be gained from simple moment equations. The two-fluid 
model (TFM) introduced by Axford Leer & Skadron (1977) and Drury & Volk (1981) 
(DV, hereafter) treats the thermal and CR populations as separate fluids coupled only 
through the hydrodynamic equations. The main effect of this coupling is a deceleration of 
the inflowing gas in front of the shock by the pressure gradient of counterstreaming CRs 
accelerated at the shock and, as a result, an enhancement of the total shock compression, 
the multiplicity of solutions and a much higher acceleration efficiency. Unfortunately, the 
underlying particle distribution implies the pressure divergence and is underdetermined in 
some other ways (see e.g., Drury 1983 ; Achtcrbcrg, Blandford fc Periwal 1984 ; Kang fc 



Jones 199~0|, and below). 



1.2. Renormalization of the two-fluid model 

A renormalization procedure to overcome the above ultraviolet divergence has been 
suggested recently by Malkov & Volk (1996) (MV96, hereafter). This theory produces 
basically the same two-fluid hydrodynamics except the renormalized CR specific heat ratio 
T instead of the usual 7 C which results from the losses. Under the assumption 7 C = 5/3 
upstream (which automatically implies that far upstream T = 5/3 as well, for the case of 
reacceleration considered in MV96), 7 C = 4/3 downstream, and in the limit pi — > oo, the 
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renormalized two-fluid model (RTFM) produces a solution qualitatively similar to that of 
the unrenormalized theory. In the case of 7 C = 4/3 upstream, when also T decreases r < j c , 
the results change dramatically. Namely, the shock compression r becomes much larger 
than the usual unrenormalized result r = 7 (for strong shocks). What happens is a very 
fast increase of the losses with decrease of 7 C which rises the compression that even tends 
to infinity (for infinitely large Mach numbers when, in addition, one takes the maximum 
possible spectral slope q — 3r/(r — 1) at the upper cut-off). This regime was not (and 
could not be) explored in MV96 since such a high compression shock requires a detailed 
information about j c , T, and R, the spectrum averaged CR-diffusivity across the shock 
transition. This would practically be equivalent to the full kinetic solution. To obtain such 
a solution was one of the main objectives of a companion paper (Malkov 1997a, paper I). 
Further motivations of the present paper will be outlined in the next subsection. 

To conclude this subsection we note that the assumption 7 C « 4/3 upstream is precisely 
what the kinetic solution obtained in paper I suggests for the case of injection in contrast to 
the case of reacceleration, 7 C = 5/3 considered in MV96. Moreover, 7 C is smaller upstream 
than downstream, again opposite to the reacceleration case. Finally, the RTFM results in 
the injection case are unacceptably sensitive to the values of 7 C and F that are not known 
to the required extent when the kinetic solution is not available. This makes the moment 
approach especially restrictive for describing namely the injection triggered acceleration 
process. As in paper I we concentrate here on this, certainly more interesting and at the 
same time more difficult case. 



1.3. 'Pathological' limits of the TFM 

There exists another difficulty of the TFM that has already been criticized in the 
literature (e.g., | Jones fc Ellison 199 1|) . Namely, the shock modification in a steady state 
occurs while particles are constantly injected at some rate, are then accelerated, and 
disappear eventually through the upper cut-off or downstream. Once the injection is 
somehow eliminated, the CR-dominated (or efficient) steady-state solution cannot be 
justified physically, and the ordinary gas shock remains the only solution possible. The 
TFM, however, still produces a CR-dominated shock that even becomes completely smooth 
beyond a certain Mach number, M > M\. In fact, it represents a fast-mode shock in a 
two-fluid hydrodynamics associated with the tenuous high pressure CR-fluid ( Ptuskin" 



1981| ). Since the number density of the CRs n c is irrelevant in the TFM, solutions that 
have a finite CR pressure (P c > 0) are formally permitted without any injection (n c = 0). 
Moreover, for an arbitrarily small nonzero injection rate there exists a critical Mach number 
M2 < 00 above which this efficient solution is the only one the TFM can offer. 

Clearly, it is difficult to judge the acceleration efficiency because the injection rate is, as 
a rule, very small and the two-fluid system, on the other hand, does not behave adequately 
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when injection vanishes. The question, however, is whether the consideration of this limit 
within the TFM is admissible. The answer is definitive not. Indeed, a fully kinetic steady 
state solution obtained in paper I revealed the following nonlinear response of the system 
to a weak injection of thermal particles. First of all, this response depends not so much on 
the injection itself but rather on the parameter A\ = u/5 = upi/po, where v and po are the 
injection rate and injection momentum, respectively, and 5 <C 1. Furthermore, the effect of 
shock modification completely disappears as Ai — > (it becomes practically insignificant 
in an abrupt manner already at v/\/~8 ^ 1). Only in the other extreme, Ai — > oo (more 
precisely, when Ai 3> M 3//4 , see also subsection Of) a unique solution that is indeed 



injection insensitive appears which is quite in the spirit of the TFM. Quantitatively, this 
solution can be very different from the respective TFM solution but for a different reason 
which is related to the subshock smoothing. The TFM is of course incapable of describing 
the dependence of its solution on Ai since it implies p\ = oo, i.e. Ai = oo already on the 
derivation level; even if v is set to zero afterwards the physically correct behavior of the 
solution with vanishing injection cannot be recovered since Ai remains infinite. Therefore 
it is useless to expect from the TFM a correct behavior at v —>■ and to criticize this model 
for the lack of it. This is beyond its validity range. 

The most dramatic consequence of this (pi — > oo) degeneracy is the subshock smoothing 
(r s = 1). The solution becomes enormously different from the kinetic solution that does not 
pass through the point r s = 1 just because of this fact. Why this is so, has been explained 
in paper I and we shall look at this problem from a different perspective in subsection |4.1| . 

Perhaps the most direct explanation why the kinetic solution differs so strongly from 
its hydrodynamic counterpart is the singular character of the underlying perturbation 
problem in the small parameter 5 = Po/pi *C 1. No matter how small it is, the efficient 
kinetic solution with 5 = (and, hence, the TFM solution) is fundamentally different from 
that with 5 > (see paper I). 

The parameters that the TFM usually operates on are the Mach number M and the 
injection rate v (or the seed particle pressure in the case of reacceleration). As we argued, 
this is not enough to describe consistently what is going on in the steady nonlinear shock 
acceleration. The RTFM introduces an additional parameter needed, pi, and accounts of 
the losses at p = p\ but then it lets p\ — > oo, v fixed. That means Ai = u/5 —* oo, and 
therefore the results are again insensitive to v when v — > since the critical information 
about Ai is lost, exactly as in the TFM. We emphasize that it is the parameter Ai that 
regulates primarily the budget of energetic particles at a shock, not v alone. That is why, 
for example, the subshock completely vanishes in the TFM as well as in the RTFM beyond 
a certain M = Mi even for v — > 0; the more important parameter A x remains infinite which 
effectively corresponds to the situation with a very strong injection. 

If we allow for time dependence on the kinetic level of description and assume a slower 
than in the Bohm limit momentum dependence of the CR diffusion coefficient, completely 
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smooth TFM stationary solutions will exist even with zero injection and with no seed 
particles upstream. This phenomenon has been explained by Drury (1983)- the high energy 
particles may be considered as those injected in the past and being then continuously 
accelerated at the shock. Unfortunately the smooth TFM solutions can tell us very little 
about a steady acceleration of CRs out of the thermal upstream plasma - if there are no 
preexisting CRs, there are no such solutions. 

Since these TFM solutions are essentially time dependent on the kinetic level (see e.g., 
[Fallc fc Giddings 1987| , [Kang fc Jones 1990| and prury, Volk fc Berezhko 1995| for relevant 



discussions), we may infer that once a natural cut-off p\ < oo exists and is reached, this 
acceleration regime will be disrupted. Indeed, these solutions correspond to the acceleration 
of CRs injected in the past whose number density virtually decreases (although being 
irrelevant in the TFM which in fact admits these pseudo-steady solutions) while the CR 
pressure remains approximately constant, being determined simply by the ram pressure of 
the inflowing gas. When particles start to leak through the upper cut-off the CR pressure 
decreases as well and the system relaxes to the ordinary gas shock which is the only steady 
state solution without injection. In general, the above TFM quasi-stationary acceleration 
scenario implies a rather low production of CRs (in terms of their number, not the energy 
density) because it operates only on initially injected particles and suppress further injection 
as soon as this solution is set up. Much more productive would be solutions that allow for 
permanent losses. This would mean, in fact, the propagation of high-energy CRs into the 
shock surroundings which decouples them from the gas flow with their replenishment due 
to the permanent injection at the subshock. But these solutions are essentially kinetic and, 
as we emphasized, fairly different from the TFM solutions. 

As it was demonstrated in paper I, given the injection rate three different solutions are 
possible for sufficiently high M and pi, Figure 1. However, only the most efficient solution 
with the highest compression ratio has been considered in detail. Accordingly, only a first 
critical injection, i.e. the injection rate v — V\ above which this solution exists along with 
the two other solutions has been calculated. The calculation of a second critical injection 
that requires an inspection of the two remaining solutions and above which the efficient 
solution is the only possible, is one of the subjects of the present paper. We consider 
the entire manifold of stationary solutions and in this context the three above-mentioned 
solutions are simply its subsets. 

In the next section we briefly review the physical formulation of the problem and 
discuss our strategy of a unified description of all the three solutions. In Sec. 3 we obtain 
both the efficient and inefficient solutions from the integral equation derived in paper I and 
consider their matching in an intermediate range. In Sec. 4 we describe the solution space 
as a whole and calculate the critical injections. We conclude this section with implications 
of its results for the TFM . Further, in Sec. 5 we speculate upon possible scenarios of time 
dependent acceleration on the basis of the emerged bifurcation picture. Sec. 6 discusses 
the results and some of their most evident consequences for calculations of the acceleration 
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efficiency in real astrophysical objects. 



2. 



Kinetic solution 



The standard formulation of the problem of diffusive shock acceleration includes the 
diffusion-convection equation for the high energy particles constrained by conservation of 
the fluxes of mass and momentum (see e.g., Prury 1983| ; |Blandford fc Eichler 1987]) . The 
physical situation that we shall consider throughout this paper is described in paper I and 
we quote only the key elements below. 



We assume that a strong CR-modified shock propagates in the positive x-direction and 
in the shock frame of reference the steady state velocity profile of the gas is specified as 
follows: U(x) = —u(x) for x > 0, u(0+) = u ; u(oo) = u\. The equation that describes 
the isotropic part of particle distribution g(x,p) at sufficiently high momenta is the 
diffusion-convection equation (Parker 1965; Gleeson & Axford 1967) which in the region 
x > takes the form 



Here k denotes the particle diffusion coefficient that is assumed to be k(p) = nop/po = K ip/Pi, 
where po and p\ are the injection and cut-off momenta, respectively and the number density 
of the CRs is normalized to Airgdp/p. In the downstream region (x < 0) we choose, as 
usual, a spatially homogeneous solution given by U{x) = —u<i = const and g(x,p) = go(p). 
Thus u 2 is the (constant) downstream speed. Eq.([L|) governs the process of stationary 
acceleration of the CRs that are drawn (injected) from the thermal plasma at p ~ p 
and leave the system at p — pi ^> p . Since the pressure of accelerated particles directly 
influences the flow profile u(x) in the shock precursor, eq.(P should be complemented by 
the conditions of the mass and momentum conservation 



2.1. Overview of the physical formulation 




(1) 



P c + pu 2 



pu 



P\U\ 

p x u\ 



(2) 
(3) 



Here p{x) is the mass density, p\ = p(oo), P c is the CR pressure 




(4) 



and no seed particles are present, i.e. P c (oo) = 0. The particle momentum p is normalized 
to mc. Eq.(Q) is written in the region x > where we have neglected the contribution 
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of the adiabatically compressed cold gas confining our consideration to sufficiently strong 
shocks with M 2 = piu\/^P g i ^> (ui/uq)' 1 , where 7 is the specific heat ratio of the plasma 
(see paper I for a detailed discussion of this approximation). 

The subshock strength can be determined from the ordinary Rankine-Hugoniot 
condition for the gas 

uo = 7 + 1 , 

s u 2 7 - 1 + 2M~ 2 W+ 1 1 

where R = Ui/uq. 



2.2. Invariant form of the solution 

It is quite clear that the formal solution to eq.(TC]) can be regularly found as a series 



in parameter {u\ — Uq)/uq if the latter is small. Fortunately (see paper I and [Malkov 
1997b| ) also in the case of a very strong modification [u\ ^> uq) the structure of the solution 
does not change significantly, provided that it is written in appropriate variables. The key 
variable here is simply the flow potential 



* = / udx (6) 
Jo 

and the approximate solution in both cases can be represented as 

g = £o(p)exp j-^-^j (7) 

If the solution is efficient (ui ^> no), (3 is given by the following relation 

1 dhigo /q\ 



3 <91np 



In the case of inefficient solution (u\ ~ u ) the /3-term in eq.(^) should be omitted 
altogether. In fact, there is no much difference between these two cases as far as the formal 
representation of the solution is concerned; as it was shown in paper I, for efficient solutions 
/3 ~ 1/6 over a broad momentum range. On the other hand the explicit dependence of these 
solutions on p and x can be vastly different in these two cases since both the flow potential 
and the downstream spectrum may differ dramatically. It is important to emphasize that 
it is this latter circumstance that constitutes the nonlinearity of the acceleration process, 
not the difference between the representations of the efficient and inefficient solutions by 
the formula (0) which is reflected in the f3 term. The calculations in paper I show that this 
P term results in a numerical factor 9 pa 1.09 in the expression for the shock compression 
ratio instead of 6 = 1, that one would get without (3 term in the solution (|7|), i.e. taking 
it in the form of inefficient solution. Thus eq.(^) may be regarded as an invariant form of 



the approximate kinetic solution. It creates a perspective of universal description of fairly 
different acceleration regimes. For this purpose we extend the technique developed in paper 
I for describing the efficient and intermediate solutions to the case of inefficient solution. 
Acting within the above-mentioned 'ten percent' accuracy (6 = 1.09 ~ 1), we will handle 
the ^-factor somewhat loosely, putting 9 = 1 in final results when we apply them to the 
inefficient part of the solution without a more accurate matching. The price that we pay 
for the universality of description of all the three branches is not high, if one realizes that 
the global bifurcation properties of the system under question are not really known even 
qualitatively. 



3. Unified description of the acceleration process 



The basis of a universal treatment of acceleration solutions is established in paper I and 
its critical step is the derivation of an integral equation for a normalized spectral function 



J(p) 



V(p) 
V(p ) 



(9) 



where V is the spectral function introduced by the following relation 



V(p) = I < 
lo- 



with 



s(p) 



K(p)V(p) 



u 2 + V(p) + -- 



3 d\np 



(10) 



(11) 



The explicit dependence of the variable s on V is not very critical here since V effectively 
cancels out in eq. ([LTD and eq.([T0|) should be regarded logically as an integral transform 
u(^f) i — > V(p) rather than an equation for V(p) given m(\E'). An independent integral 
equation for V(p) (or J) was, in turn, derived from the Bernoulli's integral (|3|). Using 
a number of approximations discussed in paper I, this equation can be written down as 
follows 



dr' 



1 



T' + T T'J{ 



T 



exp 



3 


f dr" 


e(r s -i)J 


e t"J(t")_ 



1. 



(12) 



We have used the notations 
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We use here v as the injection rate that differs from the injection rate r/ used in paper I and 
is specified as follows 

p r s mn c c 2 

v=vm- 7 — 2~ ( 15 ) 

where 

n c = 47T / g Q (p)dp/p (16) 

Jpo 

In order to close the system formed by eqs.(|5]) and ( |i~2"|) we need another equation to relate 
the three variables (z/, R, r s ) among whose only one, say v we consider as given. As an 
intermediate step we relate the precursor compression R and the spectral function J by 
inverting eq . (|T0|) : 

-— = - — r / e s *V{s)ds + Au5{V) (17) 



where 5 is a delta function corresponding to the jump of u at the subshock. Integrating 
then eq.([l7]) over \l/ between \1/ = 0+ and oo we get 



ui-u = — — J r + iO) - J r - z'O (18) 

Z7TZ J—l/e T 

We have put V"(po) ~ (see paper I for details concerning this approximation). The 
integral around the cut (— 1/e, — e) may be evaluated with the help of eq. (|l2D and the last 
equation rewrites 

t~~j = C ~U (19) 
1 - r s 1 £ 

where 

U = [ 1,£ -^—e^ (20) 



r 



= --- - — r and 0(r) = / . . (21) 



with 

() = , 

B{r s -l) ^ v ' Je t'J{t' 

In general, the system formed by eqs.(|5]), (P^D, and (|T§| ) may have multiple solutions. On 
the other hand for sufficiently small injection rates v there must always be a solution to 
eq.(|T2"D that corresponds to the test particle acceleration regime in which J — > 1 as v — ► 0. 
This solution can be written down in terms of a Neuman series in v or ( as follows 

j = i + JLi n ll±I + o(C 2 ) (22) 
£ + r 

where we have put 6 = 1 (see the preceding section) and r s = 4 for simplicity. Solution 
(|2|) is essentially perturbative and cannot describe two remaining solutions of eq.(|I"2D that 
appear beyond some v = V\ > 0. As we shall see all the three solutions may be conveniently 
described by the single valued function v = v{R) in the (R, u) plane (Figure 1). We term 
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the solution with R > R% efficient, and that with R 2 < R < R± - intermediate. It merges 
with the inefficient solution at the point v = u 2 , R = R 2 . For a fixed v € (vi, u 2 ) all the 
three solutions have different values of R and hence different subshock compression ratio r s . 

For the further bifurcation analysis the following relation appears to be useful. It may 
be derived from eq.(|l^) by multiplying both sides by (l/r 2 J) exp(Q(j)) and integrating then 
the result between e and e -1 : 

— [7 2 + (l-ft)[/ + £e^°-- = (23) 

where 

0o = ^(e- 1 ) 



Summarizing this subsection we note that eqs.fl5|), ([12]), and (|19|) form a closed system for 
describing all the three branches of the solution. 



3.1. Branches 2 and 3 

As we mentioned the branches 2 and 3 of the function -R(z/) in Figure 1 were studied 
in paper I with the emphasis on the branch 3, that corresponds to the efficient solution. In 
this approximation J ^> 1 and hence one may neglect the second term on the right hand 



side (r.h.s.) of eq.([12|) as well as the term in the exponent since ( > e. After the formal 
solution of the resulting equation for J is obtained, eq.(|T3"D for ( serves as a nonlinear 
'dispersion relation' for this system. There is a pair of solutions in the region v > v\ and 
there are no solutions if v < v\ within the approximation ( > e. Note that the inefficient 
solution (branch one, ( e) still exists and remains the only possible for v < v\. The 
above approximation becomes better and better with increasing v > v\ for the branch 3 (( 
grows with v) and it becomes worse for the branch 2, since in the latter case ( decreases 
and the neglected terms become more and more important. We obtain here the part 2-3 
of the bifurcation curve R(v) in Figure 1 direct from eq.(p3"D rather than from the explicit 
solution of eq. (|12"D as it was done in paper I. The former procedure is more convenient 
for matching with the rest of this curve. First, neglecting the second and the third terms 
eq.(|23"D simplifies to 

U * ^ (24) 

The last approximation, being strictly valid for the branch 3, is also good on some part of 
the branch 2 provided that (see eqs.(|]l 



\/C 

— a|i-n| (25) 



e 
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Inequality (^6|) may be obtained estimating the integral in eq.(pT3[) as y 2/C (paper I). These 
conditions are, generally speaking, opposite to the validity condition of the low £ expansion 
(fZ2]) (branch 1) that is 

y/C < e (27) 

This makes the analytic description of the intermediate branch 2 (Figure 1), where a/C may 
approach e, more difficult as it was explained earlier. Fortunately, the factor |1 — Q\ in 
eq.(|2"oD is rather small numerically, unless r s « 1, and there exists an intersection between 
the half-intervals ( p5j) and (|27|) . The smallness of 1 1 — Q | is ensured by the facts that 9 — 1 
is numerically small and r s is close to 4 on the branch 2 since R is appreciably smaller there 
than on the branch 3 for the same M and v. The inequality (^) is formally much more 
restrictive regarding the matching with the branch 1. Nevertheless, as we shall see in the 
next subsection, a quite accurate matching is still possible because the branch 1 can in fact 
be extended to the region (p6|). Besides that the third term in eq.(^), whose neglecting 
requires the constraint (|26|), will not significantly exceed, in fact, the last term even for 
R -C R\ (see eq.(^) below). Coming back to eq.(^) for the flow compression R we obtain 




R-l = 1 it (2£ 

We have used eq.(||) with 7 = 5/3 and eq. (|24l) . Now we are able to obtain the 2-3 piece 
of the curve v = u(R). On the branch 3 and partly on the branch 2 not far from the 



critical point (vi,Ri) we may replace the integral in eq. (|I~3|) by y2/C Eq. (|I~3|) can thus be 
transformed to the required form 

^M)4exp^- i ^ 8/3M . 2 j (29) 

where the function ((R, M) should be taken from eq.([28|). The function v is plotted against 
R for M = oo and 5 = Po/pi = 10~ 4 in Figure 2 with the dashed line. As we mentioned this 
solution fails to work for R sufficiently smaller than Ri (Figure 1) and should be replaced 
by a different solution that must be valid down to R = 1, v = 0. The overlapping region 
to match both asymptotic solutions is presumably within the interval Ri < R < R\- For 
this purpose we obtain in the next subsection an approximate solution that describes the 
branch 1 reasonably well, in particular it yields a qualitatively correct behavior for v — > 0, 
is also valid for the branch 2 and may thus be matched with the solution (f2~9|). 



3.2. Branches 1 and 2 

The key step here is to rewrite eq.(|l2|) as follows 

C /* ^ 



JU) = 1 / / * \ + 1 (30) 

KYJ eJo r(<f/) + t{4>) v ; 
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where the function r(<j>) is defined by eg. (|21|) . Our nearest goal is to work out a usable 
approximation to the solution of this equation that is good for the region of sufficiently 
small R < R\. As we have seen in the previous subsection, the branch 3 is characterized 
by the condition Q<po ~ 1 for R ~ Ri and Q<po <C 1 for R ^> R±, which facilitates the 
study of efficient solution far from the critical point R—R\. Therefore it is natural to use 
the approximation Q<f) ^> 1 for R < Ri and to examine the possibility of matching the 
obtained asymptotic result with the solution fl2"9"|). In this approximation (f20o ^ 1) the 
integral in eq.(|30|) is dominated by its upper limit and we may write it as 

J M = ST+TF + 1 (31) 

We will continue the discussion of the accuracy of this approximation for moderate values 
of Q(J)q in the next section arguing in terms of a very good matching of this solution at 
Q<po ~ 1 with the solution given in the preceding subsection for the 2-3 piece of the curve 
u(R). Substituting now the asymptotic result fl3~T|) in eq. fl21~|) we get 



where 

e = vR (33) 
Next, since the flow compression R is not very large for 1-2 branches we can simplify the 
algebra by using the approximation of infinite Mach number for these branches, i.e we put 
r s = 4 in eq.(^). This is a good approximation for R <C M 3 / 4 . Neglecting also 6 — 1 « 0.1 
compared to R - 1 £ 1, in eq.fl^), eqs.(|13p, (|IJ) and fl2|) yield 

(R - l) 2 = |e (^e-^ - l) (34) 

The last approximation is formally inaccurate for R— 1 <^ d — las the inspection of eq. (P5|) 
may show. On the other hand the qualitative behavior of v(R) remains unchanged in 
this region and is essentially v ~ R — 1. As we mentioned earlier the value 9 w 1.09 was 
calculated in paper I for the efficient (intermediate) solution whereas for the inefficient 
solution in the region R — 1 <C 1 the value = 1 is, in fact, a better choice. In this case 
eq.(|34]) is virtually also valid for small R — l provided that r s m 4, which is definitely the 
case for such a small R. Substituting O from eq. (|32"D and taking the relation Q ^ 1/8 into 
account we finally obtain 



-1 £ z i+e S _ 1 

2 + ed 



(35) 



Together with the relation ( p3[) this equation determines the 1-2 part of the curve v = u(R). 
It is plotted in Figure 2 with the solid line. The both solutions may be matched smoothly 
by the construction of an appropriate intermediate expansion. It is in fact not needed for 
our purposes since the deviation of these two solutions is numerically insignificant in an 
extended region between the intersection points. 
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4. The overall picture of stationary acceleration 

In the preceding section we have demonstrated that the method of integral equation 
developed in paper I allows one to describe the process of steady shock acceleration on a 
universal basis, which gives to the notions of efficient and inefficient acceleration a concrete 
mathematical content in the form of a relevant bifurcation analysis. Although different 
asymptotic approaches for describing efficient and inefficient regimes are applied in the 
present study, the both solutions match in the region of the intermediate acceleration 
regime very well. This is in fact an overlapping region for the two solutions in the sense 
of the obtained asymptotic results. Generally, the integral equation fll2|) describes all the 
three regimes on the same basis so that the matching procedure can be improved and a 
uniformly valid asymptotic expansion can be constructed. However the achieved accuracy 
is sufficient for our purposes. 

In fact we are already in a position to describe the stationary solution space in terms 
of system parameters and dependent variables that characterize the acceleration process. 
In principle, our parameter space is two-dimensional and contains the Mach number M 
and the cut-off momentum p\ (we will use the parameter p±/po = 5 _1 3> 1). It should 
be reminded that our consideration is restricted to the region M 3> R 5 ^ 6 , p\ ^> 1 (Sec.Q). 
A convenient dependent variable is the flow compression R that obviously signifies the 
efficiency of acceleration. In the present study, however, we add to this parameter space also 
the injection rate u, since this latter, even though being in principle calculable, may vary 
depending on the model of the subshock dissipation used (see the corresponding discussion 
in Appendix A of paper I). Furthermore, as we have seen, the solution may be conveniently 
represented in the form of a single valued function v{R) that is shown in Figure 2. Thus 
we perform our bifurcation analysis here in three dimensional parameter space. In a future 
work we will reduce the parameter space to its natural dimensionality (two) by specifying 
the injection rate given the subshock conditions. 

In fact the character of bifurcations may be seen from the surface plots v = v{R, M) at 
fixed p\ and v = v{R,pi) at fixed M, respectively. These surfaces are shown in Figure 3a,b. 
One sees that the multiplicity of the solution R = R(v) is always present for sufficiently 
large values of M and p±. We have not plotted, however, the region of lower Mach numbers 
in Figure 3a where the solution indeed becomes unique for all v. The approximation 
M — > oo adopted for describing 1-2 branches is inaccurate there and the usage of eq. (|23"D 
is required for plotting the surface v(R, M) instead of a simplified version of it given by 
eq. (|24D . This is clearly indicated by a less accurate matching of the 1-2 and 2-3 asymptotic 
solutions for smaller M in Figure 3a. The tendency to the uniqueness for lower M is 
quite obvious already from this figure. A more accurate investigation of this region can be 



conveniently done by a numerical solution of eq . ([12]) , or better, of a more general equation 
in paper I which we leave for a separate study. 



As it is seen from Figure 3a, the effect of the finite Mach number M results in the fast 
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increase of u(R, M) in the region R <; M 3//4 , see also Figure 4. Therefore, for sufficiently 
large v = const where M 3 / 4 < u/5, R scales as R ~ M 3//4 . In the opposite case, M 3 / 4 > v/5 
the compression ratio, represented as a function of M saturates at the level R ~ u/S. In 
this regime the behavior of the solution is close to that for M = oo. 



4.1. Critical injections. The TFM paradox revisited 

The next issue we address here is the calculation of critical injections V\ and v^- From 
eqs.(^) and (|29|), restricting ourselves for simplicity to the case M — > oo we calculate^] 





Ul = 2e\—; R 1 = e~ s/ \- (36) 



To calculate Vi we simplify eq.(p5|) as follows. First, as discussed in Sec. 2. 2 we may set 
6 = 1. We also assume f < 1 but 5^ > 1. Then, from eq.flSD we have (R - l) 2 = §£5~ ? . 
Substituting £ = vR for the second critical injection v% and for R2 we find the following 
results that we write with a 'logarithmic accuracy' 




2y/2 yln In | ,, , 

^2= 3 V ln i i ^2 = ^yinln 7 (37) 

Clearly, the nonlinearity of the acceleration process will be important when the CR pressure 
in front of the shock is comparable with the ram pressure. Therefore, an estimate for Vi 
can be obtained directly from the condition P c ~ p\u\. On putting g ~ Po/p, ( r s = 4) 
into P c and using the relativistic form of the partial pressure for all p in eq.([|), we 
get u 2 — l/ln(l/5). A similar approach to estimating the backreaction of injected and 
subsequently accelerated particles onto the gas flow has been also applied by previous 
authors (see e.g., Prury, Markiewicz fc Volk 1989| ). The both critical injections (eqs. (|36|, 
|37D ) are plotted versus 5 in Figure 5. One sees from Figs.l and 5 that the limit v — ► 
automatically means R — > 1 (ordinary gas shock) for all 5 > 0. Since v% vanishes with 5 
much slower than u% there is always a significant gap between these two critical injections, 
i.e. an extended region of multiple solutions even for very small values of S. 

It is interesting to trace the deformation of the bifurcation curve v(R) under the 
transition to the two-fluid limit. Strictly speaking there is no kinetic prototype for the 
smooth two- fluid efficient solution (r s = 1) as it was argued earlier. At the same time this 
is not linked directly with the fact that R fails to tend to unity as v — > in the TFM 



x The critical injection V\ obtained here without an explicit usage of the solution to eq.(p^) 
differs from the correspondent result (eq.(70)) in paper I by a numerical factor ~ 1. This 
difference was explained in detail in paper I (see eq.(57) and the text below it). 
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and we explain this last anomaly of this model first and shall return to the issue of the 
subshock smoothing later. To simplify the algebra we restrict our consideration to the case 
M = oo that is quite representative in this regard. For instance, one gets R = 7 from the 
TFM in this case, also for v = 0. Let us perform the transition 5 — ► 0+ (fluid limit) in 
the kinetic description. According to eqs.(|36|, p7|) we have v<i — > v\ — ► as 5 — ► 0. That 
means that v(R) — > uniformly in the segment The fact that R\ itself tends to 

infinity results simply from our assumption M = oo and is not important here (so does R^ 
although extremely slowly). We know from eqs.(28, |29|) that R, R\ < M 3//4 for all v > 



and R -> M 3 / 4 as v — > oo, Figure 3a. Coming back to the 5 — > limit we see that the 
function v(R, 5) — > u (R). Furthermore, iso(R) = for R G [1, and rises sharply at some 
Ro e [i?i,M 3 / 4 ). That means that if we now let the actual injection rate tend to zero from 
above, R will tend to Ro > R\ instead of R = 1 as is always the case in kinetic description 
for all 5 > 0. As we have seen Rq may be rather large whereas in the TFM Rq = 7 which is 
only because of the complete subshock smoothing. Thus, there is nothing surprising that 
R — > Rq 1 as v —>■ in the two-fluid description. From the viewpoint of the superior 
kinetic solution this transition reads 5 — > 0, v — > and yields R — ► Rq(M) ^> 1 for M ^> 1. 
The fact that i?o turns out to be equal 7 in the TFM is not important in the present 
context, as we mentioned. On inverting the sequence, i.e. letting v — > 0, 5 — * we clearly 
get R — > 1, but this is definitely beyond the TFM since the transition 5 — > has already 
been made during the derivation of the two-fluid equations. The sensitivity to the order of 
the limit transitions is a signature of the singular dependence of the kinetic solution upon 5 
at 5 = and has, at least formally, nothing to do with the TFM. 

For M < oo the above consideration needs some modification with an involved algebra 
but the main conclusion remains: since in the limit process 5 — > 0, u 2 > V\ for sufficiently 
large M, a correct behavior of R{v) at vanishing u, i.e. the transition R —>■ 1 as v — > is 
impossible within the TFM. 

Let us look now into the issue of the subshock smoothing. It was shown in paper I that 
the subshock may indeed be reduced significantly and its strength for sufficiently large M 
depends on only one parameter. This is a most important system parameter governing the 
nonlinear shock acceleration. It equals 

A = = J^M^ (38) 

The subshock strength, when it is weak, scales as (see eq.(83) in paper I) 

r s - 1 ~ 1/A, A > 1 (39) 

Since the TFM implies that S = 0, there is again nothing strange that it produces a 
completely smooth solution r s = 1. 

Therefore the TFM solution is completely understandable and in certain respects may 
serve as a limit of the full kinetic solution. The problem is that it is a singular limit and it 
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might thus formally be applied - provided that the closure parameters are found - only in a 
small (in fact measure zero) part of parameter space, namely where A -1 = 5 ■ M 3//4 /i/ = 0. 
The kinetic solution differs dramatically from the TFM solution even for 1 <C A < oo. This 
is perhaps the price we have to pay for describing an intrinsically kinetic problem quasi 
hydrodynamically. 

Turning now to the situation of finite A, i.e., p\ < oo which, as we just have seen, 
should be addressed kinetically, we note that the parameter pi may vary not only from 
one steadily accelerating shock to another but it can grow as a maximum momentum in a 
time dependent acceleration process and can thus pass through the critical values of the 
corresponding bifurcation diagram. If this time evolution is sufficiently slow the steady 
state solution with slowly varying parameters should be a good approximation provided 
that the temporal behavior of these parameters is determined. We continue our bifurcation 
analysis from this perspective in the next section. 



5. Possible scenario of time dependent acceleration 

It is convenient to argue in terms of a familiar projection method (see e.g., |Joseph 1976| 
for its hydrodynamic version). We start from a general evolution equation in a symbolic 
form that governs the acceleration process 

^ = W(w,fi) (40) 

where w is a suitable vector of state of the system that contains both the particle distribution 
g(x,p,t) and the flow structure u(x,t). Eq.(BO) implies the diffusion-convection equation 



along with the continuity and Navier-Stokes equations in appropriate representations. We 
denoted as fi the set of relevant parameters like M and v. As one of the projections of w 



(or u) we may choose R = U\/uq and the correspondent component of eq.([40"D takes the 
form of an ordinary differential equation (ODE) 

dR 

— = W 1 (R,r a , Pl ,...,fi) (41) 

Here we regard p\ as a dependent variable like R and r s and others denoted by the ellipsis 
rather than a parameter that belongs to the set fi which would be appropriate for a problem 
with the fixed cut-off momentum considered in the previous sections. One may also write 
the similar equations for r s , p\ etc. 

(IT 

= W 2 (R,r s , Pl ,...,ti) (42) 
^ = W 3 (R,r s ,p 1 ,...,f Jl ) (43) 
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What we have described in the preceding sections corresponds to W3 = as an external 
constraint and we have found in fact our steady state solutions as solutions of the system 

W 1 = 0, W 2 = (44) 

The last equation being resolved for r s simply yields the R-H relation ([|), whereas the first 
equation was studied in Sec. 3 and for a fixed v it provides the solution R = R(pi,M), 
represented e.g., in Figs. 3-4. Now it is convenient to depict this solution for different values 
of v and for a fixed M as shown in Figure 6. For p > pi cr and depending on v and p\ either 
one or three solutions occur. Let us fix v = c for which the three solutions coexist in the 
interval p± < p\ < p^\ whereas for p\ < p± (pi > Pi ) only the inefficient (efficient) 
solution occurs. It is obvious that the time scale in eq.(f4"2"|) is much shorter than that in 
eqs.(41) and (|43|) since R and p\ variations are followed by the subshock strength r s almost 



instantaneously. The latter may be then found from the equation W 2 (R } r s ,pi) = as a 
function of R and p\ and be then substituted into eqs.(f41~D and (|43|). The time scales in 



eqs.([Il]) and fl4*3| ) are generally not well separated. Nevertheless, since significant variations 
in R can take place under constant p\ (not vice versa, however), we may assume tentatively 
that pi evolves slowly compared to the transition time between equilibria of eq.(^) that 
are given by W\(R, r s (R),pi) = and plotted schematically in Figure 6 as R{p\). 

Relaxing the condition W% = suggests the following scenario of the time dependent 

acceleration. Suppose u — c — const and the cut-off momentum pi(t) increases in time 

beginning from a value that is in the interval p\ CT < p\ < P\\v) (Figure 6). As long as 

the lower branch R\(pi) is stable, the system will advance along this branch towards higher 

(2) 

Pi. When the second critical point p\ (z/) is reached, the flow profile together with the 
underlying particle distribution will restructure in such a way that a much more efficient 
acceleration regime corresponding to the upper branch R u (pi) emerges. If R u (p? } ) < M 3 / 4 
the subshock will not be reduced significantly which is the necessary condition for keeping 
injection at the same level, and the acceleration process may continue towards even higher 
Pi along the efficient branch R u . It is interesting to note that if, after this transition pi 
starts to decrease instead, there will be a hysteretic behavior in the acceleration process as 
it is shown in Figure 6. The decrease of pi may be caused by an enhance of losses due to 
e.g., parameter variations on the shock path. 

If Rn(p?) > M 3 / 4 the transition to the upper branch will be accompanied by a 
strong subshock reduction and, as a result, by a considerable variation in v. It should 
be noted, however, that if the subshock is still not very weak, say r s ^ 2, then v not 
necessarily decreases with decreasing r s (see [Malkov fc Volk 1995|) . This is mostly an effect 



of normalization of injection rate z/, the number of injected particles does decrease since 
the spectrum becomes steeper at p ~ p (see eq.(|l5D). Besides that, this consideration is 
based on a purely kinematic treatment of the leakage from the downstream medium and, 
in addition, is limited to relatively strong subshocks. When r s drops to r s ^ 1 the injection 
rate v may also become very low. This means that the solution must move to a lower level 
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v = const to remain quasi-steady. Such a process will depend very much on the injection 
model in use and on the underlying subshock dissipation mechanism which are beyond 
the scope of the present paper (see |Malkov 19974 |Malkov fc Volk 1997| ). One may only 
hypothesize that if the subshock reduction is very strong and v substantially decreases 
during this transition, the solution should either evolve being close to the critical injection 
rate v\ or, if v drops essentially below V\ the solution may return to the lower branch again 
to display thus essentially time dependent behavior. 

The above consideration implies the stability of the branch R\ so that the stability of 



R u and instability of Ri follow just from the continuity of Wi(R) in eq.fl4l|). A different 
situation occurs when the branch Ri is stable, at least for certain values of pi, and therefore 
Ri and R n are both unstable for the same p%. Then, as pi grows, the solution evolves along 

(2) 

Ri up to the point p\ — p\ , beyond which the acceleration process again cannot continue 
in a quasi steady manner. In a similar situation occurring in the TFM it is usually argued 
that the intermediate solution is unstable (DV). This is indeed typical for systems with 
S-type response curves like ours. We do not consider this question here in any detail but we 
note that an extended system of ODEs (eqs.([y], [|!| |43|) ) or even its further extensions will 



almost certainly possess a much more rich invariant manifold than the fixed points Ri : Ri 
and R u of eq.pl]). 



Interestingly, on the branches Ri and R u , R increases with M whereas it decreases with 
M on the branch R { (see e.g., Figure 4). Similarly, under fixed M and pi, Ri decreases with 
u, again opposite to Ri and R u (Figures 2, 3b, 6). This region of an anomalous behavior 
(like negative differential conductivity) may even be identifiable in steady state numerical 
simulations. It should be noted, however, that for very large Mach numbers M > i? 4//3 the 
compression ratio becomes practically independent of M on the all three branches (see also 
paper I). 



6. Discussion and conclusions 



The early works on the nonlinear ly modified CR shocks inspired the hope that, because 
of a very high acceleration efficiency in the nonlinear regime, the overall CR production 
should not be very sensitive to injection and just a qualitative understanding of this 
complicated process suffices for quantitative calculations of acceleration efficiency The 
CR dominated (efficient) solutions of the two-fluid model (DV) strongly supported this 
idea. In fact, such an optimism rests on the limited amount of energy available in the 
gas flow to be converted into CRs. Indeed, if we consider a strong shock (M ^> 1) we 
may calculate the acceleration efficiency, or the coefficient of the flow energy conversion, 
as £ con v = P c (0)/piu1 — 1 — 1/R, 6q.(0). Whenever the acceleration process is known to 
be in the nonlinear regime {R 3> 1), almost all the flow energy goes into CRs, practically 
independent of anything at all. 
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The main issue now is, under which circumstances the system may indeed be in a 
highly nonlinear acceleration regime. As we have seen, the answer to this question depends 
critically on the injection rate: infinitesimal variations of v in the vicinity of V\ or v 2 can 
result in finite (and typically very large) variations of R due to transitions between different 
branches. Also solutions belonging to the same branch are typically very sensitive to the 
injection rate v. This may be easily understood from the inspection of e.g., Figure 2 and 
Figure 3a where only on the inefficient branch R varies relatively slowly with u, while 
in the cases of intermediate and especially efficient solution R changes very rapidly with 
injection. There is of course an injection insensitive region belonging to efficient solution 
where the compression R(v) approaches its upper bound R ~ M 3 / 4 at given M (see Figure 
3a, where this region may be identified as a very sharp growth of the function v(R, M) 
in the farthermost corner of the plot). Such a behavior is caused by the requirement of 
a finite subshock. It should be noted that this scaling has indeed been observed in some 
numerical works with fixed injection rates (e.g., [Berezhko, Ksenofontov fc Yelshin 1995| ). 



However, this situation occurs at sufficiently large values of v and diminished r s and it 
is doubtful that such a high injection rate is possible at a weak subshock. For smaller u, 
when v <^ 6 ■ M 3 / 4 (the system parameter A 1) an important signature of the stationary 
acceleration is that the compression ratio is practically independent of M (see paper I for 
further details). 

Physically, the injection rate should be calculated selfconsistently using the solution 
of injection problem given subshock parameters. The solution of this problem provides 
a function v = u s (R,M). Then isolated solutions for R = R(M) might be obtained as 
intersection points of the curves v s and v (as shown in Figure 1). This solution may or 
may not be multiple depending on the character of the function u s . In any case, the 
bifurcation diagrams alone do not suffice for determining the actual acceleration efficiency. 
The calculation of the injection u s as a function of shock characteristics is equally important 
for this purpose. It should also be born in mind that the model considered here and in 
paper I gives an upper bounds to the actual acceleration efficiency. A number of not 
included factors may significantly decrease compression ratio R, acceleration efficiency and 
the spectrum hardness (see paper I for a relevant discussion) . 

Turning to the time dependent acceleration we note that unless the scenario suggested 
in the previous section is totally unrealistic, a critical quantity that would determine the 

(2) 

CR production is the cut-off momentum p\ beyond which the system jumps to the efficient 
acceleration regime. According to eq.(^) p± oc exp(l/z/). Since in the Bohm limit we may 
write Pi(t) oc t, the corresponding critical time t crit oc exp(l/z/). In an accelerating object of 
a finite life time r (e.g., supernova remnant, SNR) the main question, of course, is whether 
the condition r > £ crit is fulfilled. This is again extremely sensitive to the injection rate. 

Theoretically, the injection rate u s must not necessarily be as high as Vi for the 
acceleration process to become efficient. The condition u s > v\ could suffice provided that 
the lower branch R[ (Figure 6) looses stability for p\ <^ pi . Since for large values of pi, 
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v 2 ^ lOz/i (see Figure 5), this may determine the outcome of the acceleration process 
completely. Another possibility to overcome the high v 2 threshold is an essentially stronger 
time dependence of the acceleration process than that discussed in the preceding section. 
Basically, the bifurcation picture of this system is quite rich and promises an interesting 
dynamics. This is the more so as governing parameters are themselves subject for a 
temporal evolution. They may change significantly during the acceleration process in a 
variety of astrophysical environments. These variations may be of a quasi-external type 
like e.g., decrease of the Mach number when the shock slows down. Equally important 
may be an intrinsic variability associated with the growth of the maximum energy or 
with the heating of the upstream plasma by the CR driven turbulence ( | Volk, DruryJ 
fc McKcnzic 1984 ). Therefore, to comprehend the acceleration dynamics we must face 



the injection problem together with the physics of subshock dissipation and treat these 
problems selfconsistently with the above bifurcation analysis. The fact that this system 
displays very much hysteresis emphasizes the necessity of this approach. 

It is important to recognize that there is a serious drawback in the way to a full 
calculation of the acceleration efficiency in concrete astrophysical shocks, like e.g., SNR 
shocks. It originates from the threshold nature of the acceleration process. Indeed, the 
flow structure changes quasi- abruptly when the critical injection i/ 2 drops below z/ s (i? 2 ) as 
Pi grows (or u s (Ri) becomes smaller than v{). Since the function v s is very sensitive to 
local subshock conditions (a local orientation of the magnetic field is perhaps the most 
obvious and very important factor here), this transition occurs first at those parts of the 
shock surface where u s reaches its maximum. This must result in 'hot spots' or 'discharge' 
zones in the shock front where the acceleration becomes efficient. Then, the flow structure 
will be essentially 3-dimensional (or at least quasi 2-dimensional) , quite complicated and 
probably unsteady. The inhomogeneity of the ambient medium (see e.g., [McKee 1982; ) may 
very well result in a similar effect. Clearly, the one-dimensional calculations, even with 
a properly determined injection rate u s (R,M) may give at best an upper bound to the 
acceleration efficiency Even if the flow remains quasi-laminar the overall efficiency will be 
reduced according to the surface density of the hot spots on the shock front. Besides that 
the losses from the hot spots into the neighboring regions of inefficient acceleration may 
significantly reduce the maximum energy. A very important consequence of this would be 
the corresponding increase of the critical injection, eq. fl3"o] ) which may drive the system 
below the threshold of the efficient acceleration. 
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Fig. 1. — The nonlinear response R of an accelerating shock to the thermal injection v 
represented in the form of a single valued function v[K). Given v G (z^i, v?) there are three 
substantially different acceleration regimes. A few possible graphs of v s (R) (see Sec.|6|) are 
also drawn with the thin lines. 

Fig. 2. — a.) Response curves calculated from eq. fl29D (branches 2-3, dashed line) and from 
eq.(|35D (branches 1-2, solid line), for M = oo and 5 = 10 -4 . b.) Blow up of the peak in 
Figure 2a. 

Fig. 3. — a.) The surface of stationary solutions u(R,M) plotted for 5 = 10 -3 . b.) The 
same as in Figure 3a but v is given as a function of R and 5 for M = oo. Both surfaces are 
trimmed at high v. 

Fig. 4. — The same as in Figure 3a but this time in the form of the contourplot v = const 
and for smaller M, not shown in Figure 3a. 

Fig. 5. — Critical injections v\ and v% versus S given by eqs. (|36|) and fl3~7f) . 

Fig. 6. — Curves of constant v drawn schematically on the basis of the surface plot shown in 
Figure 3b. The bifurcation of acceleration occurs when v crosses its critical value, v cx (the 
solid curve marked by v = u CT ). For v < v cx and pi > pi cr three different acceleration regimes 
emerge for pi being within a certain interval (the solid curve marked by v = u c ). These 
regimes, Ri, Ri, R u {pi) correspond to the 1,2,3 branches in Figure 1. The two remaining 
light curves correspond to v- values c < v < u CT and v < c. 



